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TECHNIQUES FOR C07JIECTING APPROXIMATE FINITE DIFFERENCE SOLUTIONS 

David Nixon 
Ames Research Center 


SUMMARY 

A method of correcting finite-difference solutions for the effect of 
truncation error or the use of an approximate basic equation is presented. 
Applications to transonic flow problems are described and examples given. 

INTRODUCTICN 


Probably the paramount objective in computational fluid dynamics is to 
reduce the computing cost while still obtaining a solution of acceptable 
accuracy. The most direct way of reducing the cost is by algorithm improve- 
ment but there are other ways of cost reduction for practical applications. 

One way is to use a coarse discretization in the numerical procedure which 
reduces the number of data points to be computed. Another method is to use 
some approximate form of the governing equations or boundary conditions that 
leads to an easier and more rapid computation. While generally representing 
all of the essential features of the complete solution these approximate 
solutions can fail to accurately capture Important finer details. However, 
if the approximate solutions could be easily corrected in some way to account 
for these deficiencies and if the correction is universal or even partially 
universal then computing costs could be reduced without sacrificing accuracy. 
It is the derivation of such corrections that is discussed in this paper. 

The particular problems addressed in the present work are concerned with 
finite difference calculations of two-dimensional transonic flow problems, 
an Important area of aerodynamics research. A common simplification in 
transonic flow calculations is the use of the approximate transonic small 
disturbance equation with thin wing boundary conditions rather than the full 
potential equation with an exact treatment of the boundary conditions. The 
error introduced by this approximation and the means of correction are dis- 
cussed in this paper. A second form of approximation discussed is the use of 
a coarse finite difference mesh which is then corrected to give a solution 
typical of a much more refined mesh. Basically the correction method devised 
for both these cases is an extrapolation from the approximate solution to the 
"exact" solution. For the transonic flow problems considered, this correction 
is obtained only for a "similar" or "nearby" solution, for example, the flow 
around a different airfoil which has all the features of the problem to be 
corrected and which has much the same pressure distribution. In this way 
the same correction can be used for all "nearby" problems, in fact a partially 
universal correction. 


BtfcauRC of the ^’eneral nature of the preesure distributions in transonic 
flow calculations in which shock waves of shock capture regions and regions 
of rapid pressure change occur a modification of the method of strained 
coordinates developed by Nixon (ref. 1) is used. In this technique the 
coordinates are strained such that all the rapidly varying parts of the solu- 
tion are constrained to the same location which improves the range of validity 
of the corrections. 

Results for the pressure distribution around several airfoils with 
corrections for both mesh size and the use of approxiiiuite equations or 
boundary conditions are presented siiowing the satisfactory application of 
the theory. The most important restriction in the choice of both corrections 
and approximate solutions is that they rnitat represent all the essential 
features of the final solution. 


BASIC TH'WRY FOR CORRECTIONS 


Consider the mathematical problem defined by 

L(«) - 0 (1) 

and the boundary condition 

B(if) ■ F(x^) on some boundary C (2) 

where L( ) is a differential operator and B( ) is an operator which may 
be differential: F(x£) is a specified function of the coordinates x^. 

It is proposed to solve the system equations (1) and (2) by a finite 
difference method. If L^( ) and B^( ) denote difference operators corre- 
sponding to L( ) and B( ) in some computational mesh characterized by A 
then the difference equations are 

- 0 (3) 

- F(x^) (4) 

where is the converged solution of the difference equation. 

For a given difference scheme 

^ + A"R(4>^) (5) 

where R((^^) depends on the essential character of the solution, that is on 
and its derivatives; the exponent n depends on the order of accuracy 
of the difference scheme. As A 0 the ditference solution should approach 
the exact solution. If it is assumed that A is sufficiently small that all 
the essential features of the solution are represented by then R((>^) 

does not change significantly with decreasing A. In this case the difference 
solution will approach the exact solution monotonically as A 0. 
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('onniJor iu>w tlu> caito where iwo difference meHheH, diameter Ized by 
Aj and A,, are tiaed where A^ < Aj. Tlien from equation (5) 

■ -[".'"" (♦a )- ''i"“,0a,)] 

I » ^ 

where R| and R. are the two error termh In tlie finite difference approxima- 
tion. It la aHHumed that the Holut ion on the A^ meah la a fairly xood 
approxinui t ion to tlie aoliition on the fine meah A-.. Tlien the quintity 

[''/'“j(«A^)- (•.!,)] 

can be aaaumed to be amall. 

A trivial relation connecting and 41 .^ is 

% * "a, * (\ ■ *A,) <’> 

From equation (b), equation (7) the correction 

\ - *A, ■ '\(*A,)- ''l"»l(*A,)] 

which la amall and aa A^, A^ >0 tlie correction is zero. 


A correction for the coarse mesh solution could be obtained from 

equation (b) if the error terms were known — an almost 

Impuasible task a priori. However, if some similar or "nearby" solution to 
la known, denoted by and if it is assumed tliat 



'a,) ♦ sh(% - 



(*») 


where t j 


is some small parameter then an approximation to equation (7) is 





It is assumed tliat aolutlons for tlie nearby problimi are known for both Aj 
and Aj mesh systems. Since both t j and { f j ( A ,”r, ) - Aj”Rj ) 1 } are 

snail quantities then a good approximate correction is given by 


\ ^ (% - ♦a,) 


( 11 ) 


Hence, if the correction for truncation error (from the Aj mesh to the A^ 
mesh) is known for the nearby problem by direct comparison then equa- 

tion (11) will correct the solution 41 ^^ to the order of accuracy implied In 
equation ( 10 ). 
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A second form of correction le when a rlnpllfled equation or boundary 
condition is used Instead of equations (1) and (2). Thus, the system solved 
is 


L(«) - 0 (12) 

with the boundary condition 

B(«) - F(xj) on C (13) 

« A 

where L( ) and B( ) are operators that give approximate forms of the 
equations (1) and (2); the boundary C may be an approximation to the exact 
boundary C. As before, these equations will be solved by finite differences 
with the difference operator 


and boundary conditions 




0 


B^(t^) ■ F(Xj) on C 


(U) 


(15) 


where is the converged solution of the discretized problem defined by 

equations (14) and (15). It is assumed in simplifications of this nature 
that the solution of the system, equations (14) and (15), is a good approxi- 
mation to the accurate solution of equations (3) and (4). Thus 

where €2 is some small quantity. 

As before, a trivial relation between and is given by 

* (♦a - i^) (17) 

and from equation (16) 

(♦^ - ' 0(C2«^) (18) 


If a solution to a nearby problem, denoted by for both the system 

equations (3) and (4) and the simplified system equations (14) and (15) is 
known then a correction to Is given by 

‘ ♦a (♦a - ^a) 

where is the solution of the nearby problem using the simplified 

operators C^( ) and B^( ). 
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If 


(*4 - *J - (»A - M - Ja)] 

whi'ro Cj Ih a Hmall quantity tiu'n thi* form.:l accuracy of equation ( 1 * 1 ) Ih 
OCcjtj). 

Hence, if the correction for the simplified problem (from l.( ) to 
L( ), etc.) Is^knovm fur the nearby problem, then equation (19) will correct 
the solution to an order of accuracy of (tjCj). 

Since equation (11) and equation (19) are linear the principle of super- 
position can be used. Thus, if a solution to the problem defined by equa- 
tions (12) and (13) is ktu>wn for some mesh and the solutions for a 

"nearby" probli*m are known for the system equations (12) and (13) on the meshes 
characterized by A| and A^ and also for the exui. t formulation, equations (3) 
and (4), then the solution for the exact problt*m on the A 2 mesh is given 
approximately by 

A • A . 

♦Aj ■ *(% - - *4,) ‘2<» 

Hence, once the correction is known for the nearby solution only the coarse 
mesh solution is required in order to obtoln an accurate estimation of the 
exact solution. 

The principles outlined in this section can be applied to othe** approxi- 
mations in addition to those conniderod above. 

The usefulness of the correction theory depends on the range of applica- 
bility of a given correction, that is, how nuiny problems can be satisfactorily 
corrected by a given nearby solution. This can only be determined by 
exper iment . 

The basis of the present theory depends on the availability of the nearby 
81 lutlons and that the changes due to the corrections are small. The validity 
of a perturbation (correction) as outlined in this section requires that the 
changes due to the correction are smooth. In solutions with discontinuities, 
for example, transonic flows with shock waves, which is the application 
considered here, the change due to the correction is not small <n the region 
traversed by the discontinuity. In order to overcome difficulties associated 
with moving discontinuities the method of strained coordinates developed by 
Nixon (refs. 1 and 2) is used. This technique is briefly described and 
extended in the next section. 


METHOD OF STRAINED DOORDINATES 


The method of correcting approximate solutions outlined in the previous 
section requires that the changes due to the correction are small. In cases 
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containing d Iscunt Inult lea which can alter location when the correction terms 
are added, the changes are not small In the region traversed by the discontin- 
uity. A means of treating the problem of perturbations in discontinuous 
transonic flow has been described by Nixon (ref. 1); an outline is given 
below. 

Briefly, the Idea Is that the problem Is reformulated In a strained 
coordinate system In which the discontinuity remains at the same location 
throughout the perturbation and hen<'c the difficulties associated with moving 
discontinuities do not arise explicitly. The required straining is then 
found as part of the solution. The basic equations In this strained coordi- 
nate system are then perturbed about some known solution to give a linear 
equation for the perturbation quantities similar to those discussed in the 
previous section. Once the solution of the linear perturbed equation is 
known, the total perturbed solution In the physical coordinates is then 
obtained. The major restriction is that the discontinuities :.iust not be lost 
or generated during the perturbation. 

The technique described above was originally developei to treat the 
discontinuities which can lnv.tlidate a perturbation analys^s. However, the 
technique can also be applied to Increase the range of application of a valid 
perturbation analysis. An example from transonic airfoil theory concerns the 
pressure distribution around an airfoil when shock waves are present; such a 
pressure distribution is sketched in figure 1. The solid and dashed lines 
denote two nearby solutions for the pressure distribution. The solution shock 
waves are captured, that is, the expected discontinuity ?.s smeared over a few 
mesh spacings and denoted by CD and C'D*. The method ot strained coordinates, 
as given in reference 1, would strain the x-coordlnate such that the midpoints 
of the shock capture regions CD and C'D* coincide. The actual details 
of the shorx ci pture region are not considered since they are in any event 
an artificial phenomena. 

It can be seen from figure 1 that in the leading edge region the rapid 
change in Che pressure distribution can cause large pressure changes for a 
small perturbation if the location of Che pressure rise shifts slightly in the 
x-direction. This large effect seriously limits the range of validity of the 
perturbation analysis sim i* all pressure changes are assumed small. A method 
of avoiding this difficulty is to strain tht coordinates such that representa- 
tive point on the AB and A'B' curves coincide. This then increases the 
range of vaJldtcy in a similar way as the treatment of the shock waves in 
reference 1. 

A further point concerns the treatment of the shock capture regions CD 
and C'D*. In the earlier applications (refs. 1 and 2) of the theory the same 
mesh and differencial equations were used for computing the pressure distri- 
bution in all examples and hence Che shock capture characteristics were 
essentially the same for all cases. For other problems, for example, cor- 
rection of a coarse mesh solution (truncation error), the shock capture 
characteristics may differ substantially and it is desirable to correct this 
behavior. Accordingly, the coordinates are strained such that both the points 
C,C' and D,D* (the extremities of the shock capture) coincide. As before. 
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rhe actual flow detallii in thv r'^glon CO are conaidered irrelevant because 
of the artificial nature of the ahock capture. 

A more Keneral statement of the above technique is as follows. 

1. If a true discontinuity is present, the coordinate straining is such 
that the location of the discontinuities coincide. 

2. If there is a shock capturinK type of phenomena, then the straining 
is such that tiu* extremities of the capture region coincide. 

3. If large gradients are present in the solution, then the coordinate 
straining is chosen so that a representative point in the region of the large 
gradient coincides. 

These conditions constitute perhaps a large number of requirements for 
the choice of straining. However, a piecewise straining is perfectly feasible 
provided the end points of the straining (which do not move) lie in regions 
of the solution for which a small perturbation analysis is valid, for example, 
in the region BC of figure 1. 

APPLICATIONS TO TWO-DIMENSIONAL TRANSONIC FLOW COMPUIATIONS 


This section is concerned with the computation of the pressure (or 
velocity) distribution around airfoils in transonic flow using either the 
full potential equation or the transonic small disturbance equation. 

If 'f(x, y) is the perturbation velocity potential in the Cartesian 
coordinates (x, y) and p is the density then the full potential equation 
is 


Ip( 1 + ^)1 + (P^vL " 0 <21) 

AX J y 

where the Cartesian velocity components u, v are given by 

u - 1 + , V - ^ (22) 

A 7 

and 

P - [l + - u2 - v2)]^'^^’‘^ (23) 

where N is the free-stream M<ich number. Both u and v are scaled with 
respect to the free-stream velocity; y is the ratio of specific heats. 
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The transonic small disturbance equation is 


(I - 




(2A) 


where q is an arbitrary parameter. 

The tangency boundary condition for equation (21) is 


♦ (x.y*) v(x,y,) 


I + ♦ (x.y,) u(x,y^) 


y^(x) 


where y ■ apecifies the geometry of the airfoil. 

The thin airfoil boundary condl'^ions for equation (24) are 

♦y(x, tO) - y^(x) 


(25) 


(26) 


where -H), -0 denote values on the upper and lower surfaces of the airfoil 
chord line, respectively. 

In the strained coordinate system the x-coordlnate is replaced by x' 
defined as 


X ■ x' + eX| (x* ) (27) 

where Xj(x') is a straining function and c is some small parameter. The 
velocity is then given (ref. 1) by 

♦jj(x,y) - ♦^°^(x',y)[l -- cx, ^(x')] ^(x',y) (28) 

X ^ 

where , y) is a known base solution. If ^^^)(x, y) denotes a second 

known solution for a value of the parameter c then 

♦ , (x'.y) • 4 y) - <>^°\x',y)[l - ex, ,(x')]| (29) 

where 

X - x' + eXj (x' ) (30) 

Combination of equations (28) and (29) then leads to 

♦„(x,y) - ♦^?^(x',y) + i r^l'^(x,y) - ♦f?^(x',y)l (31) 

X X C L X X ^ 

where the coordinates x', x are defined by equations (27) and (30). 
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For a correction for truncation error c and c are identical aince it 
is to be the m tne change of mesh that is considered. The same result also 
applies for a correction between the small disturbance and full potential 
eouations. Also as described in section 2, the correction terms are evaluated 
from a nearby solution, for example, a different airfoil but with a similar 
velocity distribution. If the solutions to the nearby problems are denoted by 
an overbar then the corrected value of u(x, y) is givt.ii by 

u(x,y) • u^®^(x',y) [u^*^(x,y) - u^°^(x',y)j (32) 

where x* Is defined by equation (27). 

If more ttian one correction is required then 

N 

u(x,y) - u^°^(x',y) ^[u{'\xj,y) - u^°^(x',y)l (33) 

i-1 

where 


and 


+ c.x, (x') 
‘ ‘i 




EiXj 


(x') 


(34) 


(35) 


where c.x, (x*) is the straining function for the ith correction and N is 
^ ‘i 

the number of corrections. 


Choice of Straining Functions 

Two forms of straining functions are given in references I and 2. If 
there is one characteristic point, for example a shock wave. In the solution 
then a suitable (ref. 1) straining function is 


eXj(x’) - 0 , 1 < x' < 0 


(36) 


where x^ is the location of the characteristic point and e6x^ Is Its 
movement between the two known solutions. 


If there are two characteristic points, x^, Xg in the solutions then 
a suitable (ref. 2) straining is 
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xM - x’)(x’ - xA) x’ci - x*)(x' - X?) 

.x,(x-) . cA-; xXHxX-Tjr ♦ ‘K x^d . x^xx; - xj) ; 0 • ^ • 


Xj(x')*0, l<x'<0 


07) 


wlwrt* c6x! .ind • ^xi arc the movi*inentM of x' and x't rcaju'ct I vcl y , between 
. . A , v* An 

the two known NOluriona. 


If there are more than two character lat Ic polntM then a plccewlai 
veralon of equations (36) and (37) can be used. 


RKSULTS 


The idcc'H discussed above are applied to probl'.*ms In Lwo-d Imenslonal 
transonic flow, Two types of correction are examined, namely, truncation 
error and the error induced by the use of the small disturbance equation 
instead of tha full potential equation. in all of the examples tlie coordinate 
straining Includes the effect of the shock wave. Including the shock capture 
correction discussed in section 3 and the rapid pressure expansion at the 
leading edge. The choice ot the nearby solution for the correction was 
determined mainly by inspection, that Is, whether the pressure distribution 
in both cases looked similar. In all examples, an interpolation procedure Is 
used to provide data at points other than those used in the computation. 

s 

In figure 2 the pressure distribution over the upper surface of an NACA 
6AA410 airfoil at *■ 0.74 and 1.3* angle of attack is shown. The final 
result Is corrected for a mesh of 99"79 from a solution obtained on a mesh of 
38»35. The nearby solution is for the same airfoil and Mach number ^•lt at 1* 
angle of attack. All results ,tre computed using the transonic small distur- 
bance theory. It is seen that the large discrepancy between coarse and fine 
mesh calculations ne.ir the leading edge is corrected fairly well by the present 
scheme as is the shock capture region. The probable reason for the error In 
the leading edge is that the nearby result does not quite represent the 
behavior of the problem under consideration. 

An important part of the present theory is that the correction need not 
be computed for the airfoil under consideration but only for an airfoil with 
a similar form of solution. Correspondingly, some examples were computed 
usinp different Airfoils for the correction. In t igure 3 the pressure dis- 
tribution around an NACA 0012 airfoil at, ■ 0.84 at zero angle of attack 
is shown for a mesh of 99*79 corrected from a computational mesh of 38*35 
using the correction obtained for the upper surface of an NACA 64A410 airfoil 
at ■ 0.74 and 1* angle of attack. Although there is some error Just ahead 
of and Just behind the shock wave, the main discrepancy in the region of the 
leading edge has been corrected satisfactorily. A typical correction for the 
nearby solution is also shown in figure 3. 
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A Ktfc'ond type of correction discuHRcJ In thle paper is for the use of the 
traniionic Hmall diaturbance equation aa opposed to tiie full potential equation. 
An example of this type of correction la given in figure 4 where the pressure 
distribution over the upper surface of a Korn airfoil is shown for K., ■ 0.78 
and i tro angle of .ittack. This result is obtained by correcting a small 
disturbance result to correspond to a full potential result. The correction 
is obtained from cinnputatlons of both small disturbance and full potential 
equations for the Korn airfoil at M,„ ■ 0.76 and zero angle of attack. The 
agreement of the corrected result with the direct solution is satisfactory 
except behind the shock wave. In this particular example there are not many 
points on the airfoil surface behind the shock wave in the small disturbance 
solution which to a great extent accounts for the lack of resolution. 

The final example, shown in figure 5, illustrates the correction between 
small disturbance and full potential equations. In this case the correction 
is computed for a different airfoil. In figure b, the pressure distribution 
around the upper surface of an NACA 64A410 airfoil at ■ 0.72 and zero 
angle of attack is shown; the correction is obtained from solutions using 
small disturbance and full potential equations for a lOZ biconvex airfoil at 
zero angle of attack at N^. * 0.828S. The corrected result agrees satis- 
factorily with the direct solution. 

In all .>t ‘he examples computed it was essential that the approximate 
are coars* .es.i solutions and the correction solutions must capture all of the 
charactei'istic features of the problem under consideration. 


CONCLUDINC REMARKS 


A technique for correcting finite difference solutions for the effect of 
truncation error and the use of approximate equations or boundary conditions 
is derived. The correction terms need only be computed from a similar or 
nearby solution. It is essential chat both the correction and the approximate 
solutions should capture all of the characteristics of the final solution. 
Since a given correction should be applicable to at least a moderate range of 
problems the present method should be very useful In practical computations 
since only the relatively Inexpensive approximate computations are required. 

The author would like to acknowledge the help of T. L. Holst of Anes 
Research Crater who provided the solutions for Che full potential equations. 
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F-niire PrLssuro d i Ht r i mil .on around a NACA 0012 airtoil 
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Figure U.- Presfsure distribution around the upper surfait c: a Kor: L.rfc 

Mo - 0.78. 3 - 0 *. 
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FULL POTENTIAL EQUATION SOLUTION 

SMALL DISTURBANCE SOLUTION 

O CORRECTED SMALL DISTURBANCE 

SOLUTION 

V NEARBY SOLUTION (10% BICONVEX 

AIRFOIL, u • 0 . M, - 0.8285) 

SMALL DISTURBANCE 



X 


Figure 5.- Pressure distribution around the upper surface of a NACA 64A410 

airfoil; M, ■ 0.72, a ■ 0®. 
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